function [z, df, d2f] = FUNC3(a,b,c,d,e,f)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This version is checked by Jun Yu on Oct 31, 2023
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% a omega0_t
% b RL_t+1
% c H_t
% d RK_t+1
% e sigma0
% f sigma1

% checked
fun = @(x,y)exp(-(x-(-0.5.*e.^2)).^2./(2*e.^2))./(e.*sqrt(2.*pi)).*exp(-(y-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi));
ymin = @(x)log(b.*c./(d.*(1+c)))-x;

z = integral2(fun,a,100,ymin,100);

% checked
funa =@(s)exp(-(a-(-0.5.*e.^2)).^2./(2*e.^2))./(e.*sqrt(2.*pi)).*exp(-(s-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi));
da=-integral(funa,log(b*c/(d*(1+c)))-a,100);

% checked
funb = @(x)1./(b).*exp(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*exp(-(x-(-0.5.*e.^2)).^2./(2*e.^2))./(e.*sqrt(2.*pi));
db= -integral(funb,a,100);

% checked
func = @(x)1./(c.*(1+c)).*exp(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*exp(-(x-(-0.5.*e.^2)).^2./(2*e.^2))./(e.*sqrt(2.*pi));
dc=-integral(func,a,100);

% checked
fund = @(x)1./(d).*exp(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*exp(-(x-(-0.5.*e.^2)).^2./(2*e.^2))./(e.*sqrt(2.*pi));
dd=integral(fund,a,100);

% checked
fune = @(x,y)(-exp(-(x+0.5.*e.^2).^2./(2*e.^2))./(e^2.*sqrt(2.*pi))+exp(-(x+0.5.*e.^2).^2./(2*e.^2))./(sqrt(2.*pi)).*(-((x+0.5.*e.^2).*e^2-(x+0.5.*e.^2).^2)./e.^4)).*exp(-(y-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi));
de = integral2(fune,a,100,ymin,100);

% checked
fung = @(x,y)exp(-(x-(-0.5.*e.^2)).^2./(2*e.^2))./(e.*sqrt(2.*pi)).*(-exp(-(y+0.5.*f.^2).^2./(2*f.^2))./(f.^2.*sqrt(2.*pi))+exp(-(y+0.5.*f.^2).^2./(2.*f.^2))./(sqrt(2.*pi)).*(-((y+0.5.*f.^2).*f.^2-(y+0.5.*f.^2).^2)./f.^4));
dg = integral2(fung,a,100,ymin,100);

df=[da db dc dd de dg];


% second order derivatives
% a omega0_t
% b RL_t+1
% c H_t
% d RK_t+1
% e sigma0
% f sigma1
    % define the probablity density function
    function [f0] = desity(x,sigma_under)
      f0 = 1/(sigma_under*sqrt(2*pi))*exp(-(x+0.5*sigma_under^2)^2/(2*sigma_under^2));
    end
    
    % call the probability density function
    f0_omega0bar = desity(a,e);
    f1_omega0bar = desity(log(b*c/(d*(1+c)))-a,f);
    
    % first order derivative of probability density function with respect
    % to sigma
    function [f0prime] = desityprime(x,sigma_under)
      f0prime = -1/(sigma_under^2*sqrt(2*pi))*exp(-(x+0.5*sigma_under^2)^2/(2*sigma_under^2))+ 1/(sqrt(2*pi))*exp(-(x+0.5*sigma_under^2)^2/(2*sigma_under^2))*(-(((x+0.5*sigma_under^2)*sigma_under^2 -(x+0.5*sigma_under^2)^2)/(sigma_under^4)));
    end
    
    % Second order derivative of probability density function
    function [f0prime2] = desityprime2(x,sigma_under)
      f0prime2 = 2/(sigma_under.^3.*sqrt(2.*pi)).*exp(-(x+0.5.*sigma_under.^2).^2./(2.*sigma_under.^2)) ...
               + 1/(sqrt(2*pi)).*exp(-(x+0.5.*sigma_under.^2).^2./(2.*sigma_under.^2)).*((((x+0.5.*sigma_under.^2).*sigma_under.^2 -(x+0.5.*sigma_under.^2).^2)/(sigma_under.^5)))...
               + 1/(sqrt(2*pi)).*exp(-(x+0.5.*sigma_under.^2).^2./(2.*sigma_under.^2)).*(-((2.*x.*sigma_under.^2+2.*sigma_under.^4-6.*sigma_under.^2.*(x+0.5.*sigma_under.^2)+4.*(x+0.5.*sigma_under.^2).^2)./(sigma_under.^5)))...
               + 1/(sqrt(2*pi)).*exp(-(x+0.5.*sigma_under.^2).^2./(2.*sigma_under.^2)).*(-(((x+0.5.*sigma_under.^2).*sigma_under.^2 -(x+0.5.*sigma_under.^2).^2)./(sigma_under.^4))).*(-(((x+0.5*sigma_under.^2).*sigma_under.^2 -(x+0.5*sigma_under.^2).^2)/(sigma_under.^3)));
    end
    f0prime_omega0bar = exp(-(a-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))*(-(a-(-0.5.*e.^2))./(e.^2));

    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For da functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % checked     
    funaa = @(x)f0prime_omega0bar.*exp(-(x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi));
    daa1=integral(funaa,log(b*c/(d*(1+c)))-a,100);
    
    % get daa %checked, a problem found on Oct 31, 2023
    daa = -f0_omega0bar.*f1_omega0bar-daa1;
    
    % get dab %checked
    dab = f0_omega0bar.*f1_omega0bar./b;
    
    % get dac %checked
    dac = f0_omega0bar.*f1_omega0bar./(c.*(1+c));
    
    % get dad %checked
    dad = f0_omega0bar.*f1_omega0bar./(-d);
    
    % get dae %checked
    f0prime_omega0bar_sigma0 = desityprime(a,e);
    funae = @(x)f0prime_omega0bar_sigma0.*exp(-(x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi));
    dae = -integral(funae,log(b*c/(d*(1+c)))-a,100);
    
    % get dag %checked
    funag = @(x)f0_omega0bar.*(-1./(f.^2.*sqrt(2.*pi)).*exp(-(x+0.5.*f.^2).^2./(2.*f.^2))+ 1./(sqrt(2.*pi)).*exp(-(x+0.5.*f.^2).^2./(2.*f.^2)).*(-(((x+0.5.*f.^2).*f.^2 -(x+0.5.*f.^2).^2)./(f.^4))));
    dag = -integral(funag,log(b*c/(d*(1+c)))-a,100);

    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For db functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % get dba 
    dba = dab; %checked
    
    % get dbb %checked
    funbb1 = @(x)(1./b.^2).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)));
    funbb2 = @(x)(1./b).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./b;
    
    dbb = integral(funbb1,a,100)-integral(funbb2,a,100);
    
    % get dbc  %checked
    funbc = @(x)(1./b).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./(c.*(1+c));
    dbc = -integral(funbc,a,100);
    
    % get dbd %checked
    funbd = @(x)(1./b).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./(-d);
  
    dbd = -integral(funbd,a,100);
       
    % get dbe %checked
    funbe = @(x)(1./b).*exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-1./(e.^2*sqrt(2.*pi))*exp(-(x+0.5.*e.^2).^2/(2.*e.^2))+ 1./(sqrt(2.*pi))*exp(-(x+0.5.*e.^2).^2./(2.*e.^2)).*(-(((x+0.5.*e.^2).*e.^2 -(x+0.5.*e.^2).^2)./(e.^4))));
    dbe = -integral(funbe,a,100);
    
    % get dbf %checked
    funbg = @(x)(1./b).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi))).*(-1./(f.^2.*sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2))+ 1./(sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2)).*(-(((log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).*f.^2 -(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2)./(f.^4))));
    dbg = -integral(funbg,a,100);
      
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For dc functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    
    % get dca 
    dca = dac; %checked
    
    % get dcb
    dcb = dbc; %checked
    
    % get dcc; %checked
    funcc1 = @(x)(1./(c.*(1+c))).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./(c.*(1+c));
    funcc2 = @(x)((1+2.*c)./(c.^2.*(1+c).^2)).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)));

    dcc = -integral(funcc1,a,100) + integral(funcc2,a,100);
    
    % get dcd %checked
    funcd = @(x)(1./(c.*(1+c))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./(-d);
    dcd = - integral(funcd,a,100);
    
    % get dce %checked
    funce = @(x)(1./(c.*(1+c))).*exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-1./(e.^2*sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2))+ 1./(sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2)).*(-(((x+0.5.*e.^2).*e.^2 -(x+0.5.*e.^2).^2)./(e.^4))));
    dce = -integral(funce,a,100);
    
    % get dcg %checked
    funcg = @(x)(1./(c.*(1+c))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi))).*(-1./(f.^2.*sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2))+ 1./(sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2)).*(-(((log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).*f.^2 -(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2)./(f.^4))));
    dcg = -integral(funcg,a,100);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For dd functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    
    % get dda
    dda = dad; %checked
    
    % get ddb
    ddb = dbd; %checked
    
    % get ddc;
    ddc = dcd; %checked
    
    % get ddd %checked
    fundd1 = @(x)(1./d).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./(-d);
    fundd2 = @(x)(1./(d.^2)).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)));

    ddd = integral(fundd1,a,100) - integral(fundd2,a,100); 
    
    % get dde %checked
    funde = @(x)(1./d).*exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-1/(e.^2.*sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2))+ 1./(sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2)).*(-(((x+0.5.*e.^2).*e.^2 -(x+0.5.*e.^2).^2)./(e.^4))));
    dde = integral(funde,a,100);
    
    % get ddg %checked
    fundg = @(x)(1./d).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(-1/(f.^2.*sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2))+ 1./(sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2)).*(-(((log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).*f.^2 -(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2)./(f.^4))));
    ddg = integral(fundg,a,100);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For de functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    
    % get dea
    dea = dae; %checked
    
    % get deb
    deb = dbe; %checked
    
    % get dec
    dec = dce; %checked
    
    % get ded
    ded = dde; %checked
    
    % get dee; %checked
    % may be problematic
    funee = @(x,y)(desityprime2(x,e)).*exp(-(y+0.5.*f.^2).^2./(2.*f.^2))./(f.*sqrt(2.*pi));
    ymin = @(x)log(b.*c./(d.*(1+c)))-x;

    dee = integral2(funee,a,100,ymin,100);

    % get deg %checked
    funeg = @(x,y)(-1./(e.^2*sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2))+ 1./(sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2)).*(-(((x+0.5.*e.^2).*e.^2 -(x+0.5.*e.^2).^2)./(e.^4)))).*(-1./(f.^2.*sqrt(2.*pi)).*exp(-(y+0.5.*f.^2).^2./(2.*f.^2))+ 1./(sqrt(2.*pi)).*exp(-(y+0.5.*f.^2).^2./(2.*f.^2)).*(-(((y+0.5.*f.^2).*f.^2 -(y+0.5.*f.^2).^2)./(f.^4))));
    deg = integral2(funeg,a,100,ymin,100);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For dg functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    
    % get dga
    dga = dag; %checked
    
    % get dgb 
    dgb = dbg; %checked
    
    % get dgc
    dgc = dcg; %checked
    
    % get dgd
    dgd = ddg; %checked
    
    % get dge
    dge = deg; %checked
    
    % get dgg %checked
    % may be problematic
    fungg = @(x,y)exp(-(x+0.5.*e.^2).^2./(2.*e.^2))./(e.*sqrt(2.*pi)).*(desityprime2(y,f));
    ymin = @(x)log(b.*c./(d.*(1+c)))-x;

    dgg = integral2(fungg,a,100,ymin,100);  
    
    % Hessian matrix    
d2f = [daa dab dac dad dae dag; 
       dba dbb dbc dbd dbe dbg;
       dca dcb dcc dcd dce dcg;
       dda ddb ddc ddd dde ddg;
       dea deb dec ded dee deg;
       dga dgb dgc dgd dge dgg];

end

